//----------------------------------------------------------------------------
// M A R I T I M E  R E S E A R C H  I N S T I T U T E  N E T H E R L A N D S
//----------------------------------------------------------------------------
// Copyright (C) 2012 - MARIN - All rights reserved - http://www.marin.nl
//----------------------------------------------------------------------------
// Program    : mermaid
// Module     : vbm
// File       : dotp_kernels.h
// Author     : M.A. de Jong
//----------------------------------------------------------------------------

#ifndef DOTP_KERNELS_H
#define	DOTP_KERNELS_H

namespace rrbsolver
{

template <class T>
__global__ void kernel_dotp1(T *odata, const T *y, const T *x, const Grid g)
{
    int cx = g.cx;
    int cy = g.cy;
    int ld = g.nx;

    int u_r1 = BORDER_WIDTH + bx * DIM_COMPUTE_BLOCK + tx;
    int v_r1 = BORDER_WIDTH + by * DIM_COMPUTE_BLOCK + ty;
    int u_r2 = BORDER_WIDTH2 + cx + bx * DIM_COMPUTE_BLOCK + tx;
    int v_r2 = BORDER_WIDTH2 + cy + by * DIM_COMPUTE_BLOCK + ty;

    int loc_r1 = ld * v_r1 + u_r1;
    int loc_r2 = ld * v_r2 + u_r2;

    int tid = Bx * ty + tx;

    __shared__ T sm[DIM_COMPUTE_BLOCK * (DIM_COMPUTE_BLOCK / DOTP_BF)];

    T sum = 0;
    
    for (int k = 0; k < DOTP_BF; ++k) {
        sum += y[loc_r1] * x[loc_r1];
        sum += y[loc_r2] * x[loc_r2];
        loc_r1 += ld * (DIM_COMPUTE_BLOCK / DOTP_BF);
        loc_r2 += ld * (DIM_COMPUTE_BLOCK / DOTP_BF);
    }

    sm[tid] = sum;
    __syncthreads();

    for (int k = (DIM_COMPUTE_BLOCK / 2) * 
                 (DIM_COMPUTE_BLOCK / DOTP_BF); k > 16; k >>= 1) {
        if (tid < k) {
            sm[tid] += sm[tid + k];
        __syncthreads();
        }
    }

    if (tid < 16)
        sm[tid] += sm[tid + 16];
    __syncthreads();
    if (tid < 8)
        sm[tid] += sm[tid + 8];
    __syncthreads();
    if (tid < 4)
        sm[tid] += sm[tid + 4];
    __syncthreads();
    if (tid < 2)
        sm[tid] += sm[tid + 2];
    __syncthreads();
    if (tid < 1)
        sm[tid] += sm[tid + 1];
    __syncthreads();
    if (tid == 0)
        odata[by * gridDim.x + bx] = sm[tid];
}

} // namespace rrbsolver

#endif // DOTP_KERNELS_H

